{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>\n",
       "@import url('http://fonts.googleapis.com/css?family=Source+Code+Pro');\n",
       "@import url('http://fonts.googleapis.com/css?family=Vollkorn');\n",
       "@import url('http://fonts.googleapis.com/css?family=Arimo');\n",
       "@import url('http://fonts.googleapis.com/css?family=Fira_sans');\n",
       "\n",
       "    div.cell{\n",
       "        width: 900px;\n",
       "        margin-left: 0% !important;\n",
       "        margin-right: auto;\n",
       "    }\n",
       "    div.text_cell code {\n",
       "        background: transparent;\n",
       "        color: #000000;\n",
       "        font-weight: 600;\n",
       "        font-size: 11pt;\n",
       "        font-style: bold;\n",
       "        font-family:  'Source Code Pro', Consolas, monocco, monospace;\n",
       "   }\n",
       "    h1 {\n",
       "        font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "\t}\n",
       "\t\n",
       "    div.input_area {\n",
       "        background: #F6F6F9;\n",
       "        border: 1px solid #586e75;\n",
       "    }\n",
       "\n",
       "    .text_cell_render h1 {\n",
       "        font-weight: 200;\n",
       "        font-size: 30pt;\n",
       "        line-height: 100%;\n",
       "        color:#c76c0c;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 1em;\n",
       "        display: block;\n",
       "        white-space: wrap;\n",
       "        text-align: left;\n",
       "    } \n",
       "    h2 {\n",
       "        font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "        text-align: left;\n",
       "    }\n",
       "    .text_cell_render h2 {\n",
       "        font-weight: 200;\n",
       "        font-size: 16pt;\n",
       "        font-style: italic;\n",
       "        line-height: 100%;\n",
       "        color:#c76c0c;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 1.5em;\n",
       "        display: block;\n",
       "        white-space: wrap;\n",
       "        text-align: left;\n",
       "    } \n",
       "    h3 {\n",
       "        font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "    }\n",
       "    .text_cell_render h3 {\n",
       "        font-weight: 200;\n",
       "        font-size: 14pt;\n",
       "        line-height: 100%;\n",
       "        color:#d77c0c;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 2em;\n",
       "        display: block;\n",
       "        white-space: wrap;\n",
       "        text-align: left;\n",
       "    }\n",
       "    h4 {\n",
       "        font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "    }\n",
       "    .text_cell_render h4 {\n",
       "        font-weight: 100;\n",
       "        font-size: 14pt;\n",
       "        color:#d77c0c;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "        white-space: nowrap;\n",
       "    }\n",
       "    h5 {\n",
       "        font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "    }\n",
       "    .text_cell_render h5 {\n",
       "        font-weight: 200;\n",
       "        font-style: normal;\n",
       "        color: #1d3b84;\n",
       "        font-size: 16pt;\n",
       "        margin-bottom: 0em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "        white-space: nowrap;\n",
       "    }\n",
       "    div.text_cell_render{\n",
       "        font-family: 'Fira sans', verdana,arial,sans-serif;\n",
       "        line-height: 125%;\n",
       "        font-size: 115%;\n",
       "        text-align:justify;\n",
       "        text-justify:inter-word;\n",
       "    }\n",
       "    div.output_subarea.output_text.output_pyout {\n",
       "        overflow-x: auto;\n",
       "        overflow-y: scroll;\n",
       "        max-height: 50000px;\n",
       "    }\n",
       "    div.output_subarea.output_stream.output_stdout.output_text {\n",
       "        overflow-x: auto;\n",
       "        overflow-y: scroll;\n",
       "        max-height: 50000px;\n",
       "    }\n",
       "    div.output_wrapper{\n",
       "        margin-top:0.2em;\n",
       "        margin-bottom:0.2em;\n",
       "}\n",
       "\n",
       "    code{\n",
       "      font-size: 70%;\n",
       "    }\n",
       "    .rendered_html code{\n",
       "    background-color: transparent;\n",
       "    }\n",
       "    ul{\n",
       "        margin: 2em;\n",
       "    }\n",
       "    ul li{\n",
       "        padding-left: 0.5em; \n",
       "        margin-bottom: 0.5em; \n",
       "        margin-top: 0.5em; \n",
       "    }\n",
       "    ul li li{\n",
       "        padding-left: 0.2em; \n",
       "        margin-bottom: 0.2em; \n",
       "        margin-top: 0.2em; \n",
       "    }\n",
       "    ol{\n",
       "        margin: 2em;\n",
       "    }\n",
       "    ol li{\n",
       "        padding-left: 0.5em; \n",
       "        margin-bottom: 0.5em; \n",
       "        margin-top: 0.5em; \n",
       "    }\n",
       "    ul li{\n",
       "        padding-left: 0.5em; \n",
       "        margin-bottom: 0.5em; \n",
       "        margin-top: 0.2em; \n",
       "    }\n",
       "    a:link{\n",
       "       font-weight: bold;\n",
       "       color:#447adb;\n",
       "    }\n",
       "    a:visited{\n",
       "       font-weight: bold;\n",
       "       color: #1d3b84;\n",
       "    }\n",
       "    a:hover{\n",
       "       font-weight: bold;\n",
       "       color: #1d3b84;\n",
       "    }\n",
       "    a:focus{\n",
       "       font-weight: bold;\n",
       "       color:#447adb;\n",
       "    }\n",
       "    a:active{\n",
       "       font-weight: bold;\n",
       "       color:#447adb;\n",
       "    }\n",
       "    .rendered_html :link {\n",
       "       text-decoration: underline; \n",
       "    }\n",
       "    .rendered_html :hover {\n",
       "       text-decoration: none; \n",
       "    }\n",
       "    .rendered_html :visited {\n",
       "      text-decoration: none;\n",
       "    }\n",
       "    .rendered_html :focus {\n",
       "      text-decoration: none;\n",
       "    }\n",
       "    .rendered_html :active {\n",
       "      text-decoration: none;\n",
       "    }\n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "    } \n",
       "    hr {\n",
       "      color: #f3f3f3;\n",
       "      background-color: #f3f3f3;\n",
       "      height: 1px;\n",
       "    }\n",
       "    blockquote{\n",
       "      display:block;\n",
       "      background: #fcfcfc;\n",
       "      border-left: 5px solid #c76c0c;\n",
       "      font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "      width:680px;\n",
       "      padding: 10px 10px 10px 10px;\n",
       "      text-align:justify;\n",
       "      text-justify:inter-word;\n",
       "      }\n",
       "      blockquote p {\n",
       "        margin-bottom: 0;\n",
       "        line-height: 125%;\n",
       "        font-size: 100%;\n",
       "      }\n",
       "</style>\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"]\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    scale:100,\n",
       "                        availableFonts: [],\n",
       "                        preferredFont:null,\n",
       "                        webFont: \"TeX\",\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# format the book\n",
    "%matplotlib inline\n",
    "import sys\n",
    "from __future__ import division, print_function\n",
    "import sys\n",
    "sys.path.insert(0,'../code')\n",
    "import book_format\n",
    "book_format.load_style('../code')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Evidence\n",
    "\n",
    "## Interpreting SAT scores\n",
    "\n",
    "Suppose you are the Dean of Admission at a small engineering college in\n",
    "Massachusetts, and you are considering two candidates, Alice and Bob,\n",
    "whose qualifications are similar in many ways, with the exception that\n",
    "Alice got a higher score on the Math portion of the SAT, a standardized\n",
    "test intended to measure preparation for college-level work in\n",
    "mathematics.\n",
    "\n",
    "If Alice got 780 and Bob got a 740 (out of a possible 800), you might\n",
    "want to know whether that difference is evidence that Alice is better\n",
    "prepared than Bob, and what the strength of that evidence is.\n",
    "\n",
    "Now in reality, both scores are very good, and both candidates are\n",
    "probably well prepared for college math. So the real Dean of Admission\n",
    "would probably suggest that we choose the candidate who best\n",
    "demonstrates the other skills and attitudes we look for in students. But\n",
    "as an example of Bayesian hypothesis testing, let’s stick with a\n",
    "narrower question: “How strong is the evidence that Alice is better\n",
    "prepared than Bob?”\n",
    "\n",
    "To answer that question, we need to make some modeling decisions. I’ll\n",
    "start with a simplification I know is wrong; then we’ll come back and\n",
    "improve the model. I pretend, temporarily, that all SAT questions are\n",
    "equally difficult. Actually, the designers of the SAT choose questions\n",
    "with a range of difficulty, because that improves the ability to measure\n",
    "statistical differences between test-takers.\n",
    "\n",
    "But if we choose a model where all questions are equally difficult, we\n",
    "can define a characteristic, `p_correct`, for each test-taker, which is\n",
    "the probability of answering any question correctly. This simplification\n",
    "makes it easy to compute the likelihood of a given score."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The scale\n",
    "\n",
    "In order to understand SAT scores, we have to understand the scoring and\n",
    "scaling process. Each test-taker gets a raw score based on the number of\n",
    "correct and incorrect questions. The raw score is converted to a scaled\n",
    "score in the range 200–800.\n",
    "\n",
    "In 2009, there were 54 questions on the math SAT. The raw score for each\n",
    "test-taker is the number of questions answered correctly minus a penalty\n",
    "of $1/4$ point for each question answered incorrectly.\n",
    "\n",
    "The College Board, which administers the SAT, publishes the map from raw\n",
    "scores to scaled scores. I have downloaded that data and wrapped it in\n",
    "an Interpolator object that provides a forward lookup (from raw score to\n",
    "scaled) and a reverse lookup (from scaled score to raw).\n",
    "\n",
    "You can download the code for this example from\n",
    "<http://thinkbayes.com/sat.py>. For more information see\n",
    "Section [download]."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The prior\n",
    "\n",
    "The College Board also publishes the distribution of scaled scores for\n",
    "all test-takers. If we convert each scaled score to a raw score, and\n",
    "divide by the number of questions, the result is an estimate of\n",
    "`p_correct`. So we can use the distribution of raw scores to model the\n",
    "prior distribution of `p_correct`.\n",
    "\n",
    "Here is the code that reads and processes the data:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```python\n",
    "class Exam(object):\n",
    "\n",
    "    def __init__(self):\n",
    "        self.scale = ReadScale()\n",
    "        scores = ReadRanks()\n",
    "        score_pmf = thinkbayes.MakePmfFromDict(dict(scores))\n",
    "        self.raw = self.ReverseScale(score_pmf)\n",
    "        self.prior = DivideValues(raw, 54)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Exam` encapsulates the information we have about the exam.\n",
    "`ReadScale` and `ReadRanks` read files and return\n",
    "objects that contain the data: `self.scale` is the\n",
    "`Interpolator` that converts from raw to scaled scores and\n",
    "back; `scores` is a list of (score, frequency) pairs.\n",
    "\n",
    "`score_pmf` is the Pmf of scaled scores. `self.raw` is the\n",
    "Pmf of raw scores, and `self.prior` is the Pmf of\n",
    "`p_correct`.\n",
    "\n",
    "![Prior distribution of `p\\_correct` for SAT\n",
    "test-takers.](figs/sat_prior.pdf)\n",
    "\n",
    "[fig.satprior]\n",
    "\n",
    "Figure [fig.satprior] shows the prior distribution of `p_correct`. This\n",
    "distribution is approximately Gaussian, but it is compressed at the\n",
    "extremes. By design, the SAT has the most power to discriminate between\n",
    "test-takers within two standard deviations of the mean, and less power\n",
    "outside that range.\n",
    "\n",
    "For each test-taker, I define a Suite called `Sat` that\n",
    "represents the distribution of `p_correct`. Here’s the definition:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```python\n",
    "class Sat(thinkbayes.Suite):\n",
    "\n",
    "    def __init__(self, exam, score):\n",
    "        thinkbayes.Suite.__init__(self)\n",
    "\n",
    "        self.exam = exam\n",
    "        self.score = score\n",
    "\n",
    "        # start with the prior distribution\n",
    "        for p_correct, prob in exam.prior.Items():\n",
    "            self.Set(p_correct, prob)\n",
    "\n",
    "        # update based on an exam score\n",
    "        self.Update(score)\n",
    "        ```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`__init__` takes an Exam object and a scaled score. It makes a copy of\n",
    "the prior distribution and then updates itself based on the exam score.\n",
    "\n",
    "As usual, we inherit `Update` from `Suite` and\n",
    "provide `Likelihood`:\n",
    "```python\n",
    "        def Likelihood(self, data, hypo):\n",
    "            p_correct = hypo\n",
    "            score = data\n",
    "\n",
    "            k = self.exam.Reverse(score)\n",
    "            n = self.exam.max_score\n",
    "            like = thinkbayes.EvalBinomialPmf(k, n, p_correct)\n",
    "            return like\n",
    "```\n",
    "\n",
    "`hypo` is a hypothetical value of `p_correct`, and\n",
    "`data` is a scaled score.\n",
    "\n",
    "To keep things simple, I interpret the raw score as the number of\n",
    "correct answers, ignoring the penalty for wrong answers. With this\n",
    "simplification, the likelihood is given by the binomial distribution,\n",
    "which computes the probability of $k$ correct responses out of $n$\n",
    "questions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Posterior\n",
    "\n",
    "![Posterior distributions of `p\\_correct` for Alice and\n",
    "Bob.](figs/sat_posteriors_p_corr.pdf)\n",
    "\n",
    "[fig.satposterior1]\n",
    "\n",
    "Figure [fig.satposterior1] shows the posterior distributions of\n",
    "`p_correct` for Alice and Bob based on their exam scores. We can see\n",
    "that they overlap, so it is possible that `p_correct` is actually higher\n",
    "for Bob, but it seems unlikely.\n",
    "\n",
    "Which brings us back to the original question, “How strong is the\n",
    "evidence that Alice is better prepared than Bob?” We can use the\n",
    "posterior distributions of `p_correct` to answer this question.\n",
    "\n",
    "To formulate the question in terms of Bayesian hypothesis testing, I\n",
    "define two hypotheses:\n",
    "\n",
    "-   $A$: `p_correct` is higher for Alice than for Bob.\n",
    "\n",
    "-   $B$: `p_correct` is higher for Bob than for Alice.\n",
    "\n",
    "To compute the likelihood of $A$, we can enumerate all pairs of values\n",
    "from the posterior distributions and add up the total probability of the\n",
    "cases where `p_correct` is higher for Alice than for Bob. And we already\n",
    "have a function, `thinkbayes.PmfProbGreater`, that does that.\n",
    "\n",
    "So we can define a Suite that computes the posterior probabilities of\n",
    "$A$ and $B$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import thinkbayes\n",
    "\n",
    "class TopLevel(thinkbayes.Suite):\n",
    "\n",
    "    def Update(self, data):\n",
    "        a_sat, b_sat = data\n",
    "\n",
    "        a_like = thinkbayes.PmfProbGreater(a_sat, b_sat)\n",
    "        b_like = thinkbayes.PmfProbLess(a_sat, b_sat)\n",
    "        c_like = thinkbayes.PmfProbEqual(a_sat, b_sat)\n",
    "\n",
    "        a_like += c_like / 2\n",
    "        b_like += c_like / 2\n",
    "\n",
    "        self.Mult('A', a_like)\n",
    "        self.Mult('B', b_like)\n",
    "\n",
    "        self.Normalize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Usually when we define a new Suite, we inherit `Update` and\n",
    "provide `Likelihood`. In this case I override\n",
    "`Update`, because it is easier to evaluate the likelihood of\n",
    "both hypotheses at the same time.\n",
    "\n",
    "The data passed to `Update` are Sat objects that represent\n",
    "the posterior distributions of `p_correct`.\n",
    "\n",
    "`a_like` is the total probability that `p_correct` is higher for Alice;\n",
    "`b_like` is that probability that it is higher for Bob.\n",
    "\n",
    "`c_like` is the probability that they are “equal,” but this equality is\n",
    "an artifact of the decision to model `p_correct` with a set of discrete\n",
    "values. If we use more values, `c_like` is smaller, and in the extreme,\n",
    "if `p_correct` is continuous, `c_like` is zero. So I treat `c_like` as a\n",
    "kind of round-off error and split it evenly between `a_like` and\n",
    "`b_like`.\n",
    "\n",
    "Here is the code that creates `TopLevel` and updates it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A 0.793123245041\n",
      "B 0.206876754959\n"
     ]
    }
   ],
   "source": [
    "from sat import Exam, Sat\n",
    "\n",
    "exam = Exam()\n",
    "a_sat = Sat(exam, 780)\n",
    "b_sat = Sat(exam, 740)\n",
    "\n",
    "top = TopLevel('AB')\n",
    "top.Update((a_sat, b_sat))\n",
    "top.Print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The likelihood of $A$ is 0.79 and the likelihood of $B$ is 0.21. The\n",
    "likelihood ratio (or Bayes factor) is 3.8, which means that these test\n",
    "scores are evidence that Alice is better than Bob at answering SAT\n",
    "questions. If we believed, before seeing the test scores, that $A$ and\n",
    "$B$ were equally likely, then after seeing the scores we should believe\n",
    "that the probability of $A$ is 79%, which means there is still a 21%\n",
    "chance that Bob is actually better prepared."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A better model\n",
    "\n",
    "Remember that the analysis we have done so far is based on the\n",
    "simplification that all SAT questions are equally difficult. In reality,\n",
    "some are easier than others, which means that the difference between\n",
    "Alice and Bob might be even smaller.\n",
    "\n",
    "But how big is the modeling error? If it is small, we conclude that the\n",
    "first model—based on the simplification that all questions are equally\n",
    "difficult—is good enough. If it’s large, we need a better model.\n",
    "\n",
    "In the next few sections, I develop a better model and discover (spoiler\n",
    "alert!) that the modeling error is small. So if you are satisfied with\n",
    "the simple mode, you can skip to the next chapter. If you want to see\n",
    "how the more realistic model works, read on...\n",
    "\n",
    "-   Assume that each test-taker has some degree of\n",
    "    `efficacy`, which measures their ability to answer SAT\n",
    "    questions.\n",
    "\n",
    "-   Assume that each question has some level of `difficulty`.\n",
    "\n",
    "-   Finally, assume that the chance that a test-taker answers a question\n",
    "    correctly is related to `efficacy` and\n",
    "    `difficulty` according to this function:\n",
    "\n",
    "        def ProbCorrect(efficacy, difficulty, a=1):\n",
    "            return 1 / (1 + math.exp(-a * (efficacy - difficulty)))\n",
    "\n",
    "This function is a simplified version of the curve used in `**item\n",
    "response theory**`, which you can read about at\n",
    "<http://en.wikipedia.org/wiki/Item_response_theory>.\n",
    "`efficacy` and `difficulty` are considered to be\n",
    "on the same scale, and the probability of getting a question right\n",
    "depends only on the difference between them.\n",
    "\n",
    "When `efficacy` and `difficulty` are equal, the\n",
    "probability of getting the question right is 50%. As\n",
    "`efficacy` increases, this probability approaches 100%. As it\n",
    "decreases (or as `difficulty` increases), the probability\n",
    "approaches 0%.\n",
    "\n",
    "Given the distribution of `efficacy` across test-takers and\n",
    "the distribution of `difficulty` across questions, we can\n",
    "compute the expected distribution of raw scores. We’ll do that in two\n",
    "steps. First, for a person with given `efficacy`, we’ll\n",
    "compute the distribution of raw scores."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def PmfCorrect(efficacy, difficulties):\n",
    "    pmf0 = thinkbayes.Pmf([0])\n",
    "\n",
    "    ps = [ProbCorrect(efficacy, diff) for diff in difficulties]\n",
    "    pmfs = [BinaryPmf(p) for p in ps]\n",
    "    dist = sum(pmfs, pmf0)\n",
    "    return dist"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`difficulties` is a list of difficulties, one for each\n",
    "question. `ps` is a list of probabilities, and\n",
    "`pmfs` is a list of two-valued Pmf objects; here’s the\n",
    "function that makes them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def BinaryPmf(p):\n",
    "    pmf = thinkbayes.Pmf()\n",
    "    pmf.Set(1, p)\n",
    "    pmf.Set(0, 1-p)\n",
    "    return pmf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`dist` is the sum of these Pmfs. Remember from\n",
    "Section [addends] that when we add up Pmf objects, the result is the\n",
    "distribution of the sums. In order to use Python’s `sum` to\n",
    "add up Pmfs, we have to provide `pmf0` which is the identity\n",
    "for Pmfs, so `pmf + pmf0` is always `pmf`.\n",
    "\n",
    "If we know a person’s efficacy, we can compute their distribution of raw\n",
    "scores. For a group of people with a different efficacies, the resulting\n",
    "distribution of raw scores is a mixture. Here’s the code that computes\n",
    "the mixture:\n",
    "\n",
    "```python\n",
    "# class Exam:\n",
    "    def MakeRawScoreDist(self, efficacies):\n",
    "        pmfs = thinkbayes.Pmf()\n",
    "        for efficacy, prob in efficacies.Items():\n",
    "            scores = PmfCorrect(efficacy, self.difficulties)\n",
    "            pmfs.Set(scores, prob)\n",
    "\n",
    "        mix = thinkbayes.MakeMixture(pmfs)\n",
    "        return mix\n",
    "```\n",
    "`MakeRawScoreDist` takes `efficacies`, which is a\n",
    "Pmf that represents the distribution of efficacy across test-takers. I\n",
    "assume it is Gaussian with mean 0 and standard deviation 1.5. This\n",
    "choice is mostly arbitrary. The probability of getting a question\n",
    "correct depends on the difference between efficacy and difficulty, so we\n",
    "can choose the units of efficacy and then calibrate the units of\n",
    "difficulty accordingly.\n",
    "\n",
    "`pmfs` is a meta-Pmf that contains one Pmf for each level of\n",
    "efficacy, and maps to the fraction of test-takers at that level.\n",
    "`MakeMixture` takes the meta-pmf and computes the\n",
    "distribution of the mixture (see Section [mixture])."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calibration\n",
    "\n",
    "If we were given the distribution of difficulty, we could use\n",
    "`MakeRawScoreDist` to compute the distribution of raw scores. But for us\n",
    "the problem is the other way around: we are given the distribution of\n",
    "raw scores and we want to infer the distribution of difficulty.\n",
    "\n",
    "![Actual distribution of raw scores and a model to fit\n",
    "it.](figs/sat_calibrate.pdf)\n",
    "\n",
    "[fig.satcalibrate]\n",
    "\n",
    "I assume that the distribution of difficulty is uniform with parameters\n",
    "`center` and `width`.\n",
    "`MakeDifficulties` makes a list of difficulties with these\n",
    "parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def MakeDifficulties(center, width, n):\n",
    "    low, high = center-width, center+width\n",
    "    return numpy.linspace(low, high, n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By trying out a few combinations, I found that `center=-0.05`\n",
    "and `width=1.8` yield a distribution of raw scores similar to\n",
    "the actual data, as shown in Figure [fig.satcalibrate].\n",
    "\n",
    "So, assuming that the distribution of difficulty is uniform, its range\n",
    "is approximately `-1.85` to `1.75`, given that\n",
    "efficacy is Gaussian with mean 0 and standard deviation 1.5.\n",
    "\n",
    "The following table shows the range of `ProbCorrect` for\n",
    "test-takers at different levels of efficacy:\n",
    "\n",
    "              Difficulty\n",
    "    Efficacy  -1.85  -0.05  1.75\n",
    "        3.00   0.99  0.95  0.78\n",
    "        1.50   0.97  0.82  0.44\n",
    "        0.00   0.86  0.51  0.15\n",
    "       -1.50   0.59  0.19  0.04\n",
    "       -3.00   0.24  0.05  0.01\n",
    "\n",
    "Someone with efficacy 3 (two standard deviations above the mean) has a\n",
    "99% chance of answering the easiest questions on the exam, and a 78%\n",
    "chance of answering the hardest. On the other end of the range, someone\n",
    "two standard deviations below the mean has only a 24% chance of\n",
    "answering the easiest questions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Posterior distribution of efficacy\n",
    "\n",
    "![Posterior distributions of efficacy for Alice and\n",
    "Bob.](figs/sat_posteriors_eff.pdf)\n",
    "\n",
    "[fig.satposterior2]\n",
    "\n",
    "Now that the model is calibrated, we can compute the posterior\n",
    "distribution of efficacy for Alice and Bob. Here is a version of the Sat\n",
    "class that uses the new model:\n",
    "\n",
    "```python\n",
    "class Sat2(thinkbayes.Suite):\n",
    "\n",
    "    def __init__(self, exam, score):\n",
    "        self.exam = exam\n",
    "        self.score = score\n",
    "\n",
    "        # start with the Gaussian prior\n",
    "        efficacies = thinkbayes.MakeGaussianPmf(0, 1.5, 3)\n",
    "        thinkbayes.Suite.__init__(self, efficacies)\n",
    "\n",
    "        # update based on an exam score\n",
    "        self.Update(score)\n",
    "```\n",
    "\n",
    "`Update` invokes `Likelihood`, which computes the likelihood of a given\n",
    "test score for a hypothetical level of efficacy.\n",
    "\n",
    "```python\n",
    "    def Likelihood(self, data, hypo):\n",
    "        efficacy = hypo\n",
    "        score = data\n",
    "        raw = self.exam.Reverse(score)\n",
    "\n",
    "        pmf = self.exam.PmfCorrect(efficacy)\n",
    "        like = pmf.Prob(raw)\n",
    "        return like\n",
    "```\n",
    "\n",
    "`pmf` is the distribution of raw scores for a test-taker with\n",
    "the given efficacy; `like` is the probability of the observed\n",
    "score.\n",
    "\n",
    "Figure [fig.satposterior2] shows the posterior distributions of efficacy\n",
    "for Alice and Bob. As expected, the location of Alice’s distribution is\n",
    "farther to the right, but again there is some overlap.\n",
    "\n",
    "Using `TopLevel` again, we compare $A$, the hypothesis that\n",
    "Alice’s efficacy is higher, and $B$, the hypothesis that Bob’s is\n",
    "higher. The likelihood ratio is 3.4, a bit smaller than what we got from\n",
    "the simple model (3.8). So this model indicates that the data are\n",
    "evidence in favor of $A$, but a little weaker than the previous\n",
    "estimate.\n",
    "\n",
    "If our prior belief is that $A$ and $B$ are equally likely, then in\n",
    "light of this evidence we would give $A$ a posterior probability of 77%,\n",
    "leaving a 23% chance that Bob’s efficacy is higher."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Predictive distribution\n",
    "\n",
    "The analysis we have done so far generates estimates for Alice and Bob’s\n",
    "efficacy, but since efficacy is not directly observable, it is hard to\n",
    "validate the results.\n",
    "\n",
    "To give the model predictive power, we can use it to answer a related\n",
    "question: “If Alice and Bob take the math SAT again, what is the chance\n",
    "that Alice will do better again?”\n",
    "\n",
    "We’ll answer this question in two steps:\n",
    "\n",
    "-   We’ll use the posterior distribution of efficacy to generate a\n",
    "    predictive distribution of raw score for each test-taker.\n",
    "\n",
    "-   We’ll compare the two predictive distributions to compute the\n",
    "    probability that Alice gets a higher score again.\n",
    "\n",
    "We already have most of the code we need. To compute the predictive\n",
    "distributions, we can use `MakeRawScoreDist` again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "exam = Exam()\n",
    "a_sat = Sat(exam, 780)\n",
    "b_sat = Sat(exam, 740)\n",
    "\n",
    "a_pred = exam.MakeRawScoreDist(a_sat)\n",
    "b_pred = exam.MakeRawScoreDist(b_sat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we can find the likelihood that Alice does better on the second\n",
    "test, Bob does better, or they tie:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "a_like = thinkbayes.PmfProbGreater(a_pred, b_pred)\n",
    "b_like = thinkbayes.PmfProbLess(a_pred, b_pred)\n",
    "c_like = thinkbayes.PmfProbEqual(a_pred, b_pred)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The probability that Alice does better on the second exam is 63%, which\n",
    "means that Bob has a 37% chance of doing as well or better.\n",
    "\n",
    "Notice that we have more confidence about Alice’s efficacy than we do\n",
    "about the outcome of the next test. The posterior odds are 3:1 that\n",
    "Alice’s efficacy is higher, but only 2:1 that Alice will do better on\n",
    "the next exam."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discussion\n",
    "\n",
    "![Joint posterior distribution of `p\\_correct` for Alice and\n",
    "Bob.](figs/sat_joint.pdf)\n",
    "\n",
    "[fig.satjoint]\n",
    "\n",
    "We started this chapter with the question, “How strong is the evidence\n",
    "that Alice is better prepared than Bob?” On the face of it, that sounds\n",
    "like we want to test two hypotheses: either Alice is more prepared or\n",
    "Bob is.\n",
    "\n",
    "But in order to compute likelihoods for these hypotheses, we have to\n",
    "solve an estimation problem. For each test-taker we have to find the\n",
    "posterior distribution of either `p_correct` or `efficacy`.\n",
    "\n",
    "Values like this are called `**nuisance parameters**` because\n",
    "we don’t care what they are, but we have to estimate them to answer the\n",
    "question we care about.\n",
    "\n",
    "One way to visualize the analysis we did in this chapter is to plot the\n",
    "space of these parameters. `thinkbayes.MakeJoint` takes two Pmfs,\n",
    "computes their joint distribution, and returns a joint pmf of each\n",
    "possible pair of values and its probability."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def MakeJoint(pmf1, pmf2):\n",
    "    joint = Joint()\n",
    "    for v1, p1 in pmf1.Items():\n",
    "        for v2, p2 in pmf2.Items():\n",
    "            joint.Set((v1, v2), p1 * p2)\n",
    "    return joint"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This function assumes that the two distributions are independent.\n",
    "\n",
    "Figure [fig.satjoint] shows the joint posterior distribution of\n",
    "`p_correct` for Alice and Bob. The diagonal line indicates the part of\n",
    "the space where `p_correct` is the same for Alice and Bob. To the right\n",
    "of this line, Alice is more prepared; to the left, Bob is more prepared.\n",
    "\n",
    "In `TopLevel.Update`, when we compute the likelihoods of $A$\n",
    "and $B$, we add up the probability mass on each side of this line. For\n",
    "the cells that fall on the line, we add up the total mass and split it\n",
    "between $A$ and $B$.\n",
    "\n",
    "The process we used in this chapter—estimating nuisance parameters in\n",
    "order to evaluate the likelihood of competing hypotheses—is a common\n",
    "Bayesian approach to problems like this."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
